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In this paper we study decaying turbulence in fixed and rotating boxes in two dimen- 
sions using the particle method SPH. The boundaries are specified by boundary force 
particles, and the turbulence is initiated by a set of gaussian vortices. In the case of 
fixed boxes we recover the results of Clercx and his colleagues obtained using both a 
high accuracy spectral method and experiments. Our results for fixed boxes are also 
in close agreement with those of Monaghan 1 and Robinson and MonagharP obtained 
using SPH. A feature of decaying turbulence in no-slip, square, fixed boundaries is 
that the angular momentum of the fluid varies with time because of the reaction 
on the fluid of the viscous stresses on the boundary. We find that when the box 
is allowed to rotate freely, so that the total angular momentum of box and fluid is 
constant, the change in the angular momentum of the fluid is a factor ~ 500 smaller 
than is the case for the fixed box, and the final vorticity distribution is different. We 
also simulate the behaviour of the turbulence when the box is forced to rotate with 
small and large Rossby number, and the turbulence is initiated by gaussian vortices 
as before. If the rotation of the box is maintained after the turbulence is initiated we 
find that in the rotating frame the decay of kinetic energy, enstrophy and the vortex 
structure is insensitive to the angular velocity of the box. On the other hand, If the 
box is allowed to rotate freely after the turbulence is initiated, the evolved vortex 
structure is completely different. 

Keywords: SPH: Turbulence, rotation, no-slip boundaries 



I. INTRODUCTION 



a ' Electronic mail: 
b ' Electronic mail: 



alireza.valizadeh@monash.edu. 



joe. monaghan@monash.edu. 



2 



Smoothed particle hydrodynamics is now widely used in computational fluid dynamics 
especially in problems involving breaking waves, and free surfaces disrupted by the impact of 
rigid bodies. In many of these systems the flows are initially laminar but develop turbulence 
as they evolve. A common example is the simulation of a laboratory dam break which 
becomes turbulent after the head of the flow hits an end wall where it forms a return wave 
that plunges into the incoming fluicP. In order to have confidence in the predicted evolution 
of the flow it is necessary to determine how accurately SPH can simulate turbulence. 

Early studies of turbulence using SPH made use of known sub-grid models applied to 
two dimensional flows. For example Shao and GotafP and Shao et alP used a 2D turbulent 
viscous stress based on the Smagorinsky model. They applied their model to the turbulence 
generated by a succession of cnoidal waves breaking on a linear ramp and used extensive 
phase averaging and other approximations (for example multiplying by a factor 4/3 to bring 
their results into better agreement with the experiments that are in 3D). They obtained 
results in reasonable agreement with experiment although SPH simulation without the tur- 
bulence model gave very similar results. Dalrymple and Rogers^ also used a Smagorinsky 
model for 2D turbulence. Violeau and IssaP studied 2D turbulence and compared the k — e, 
the EARSM (explicit algebraic Reynolds stress model), and an SPH version of Large Eddy 
Simulation. The comparison was made for a dam break, but dam breaks have the disad- 
vantage that while they are initially nearly 2D, they become 3D when the fluid strikes the 
wall of the experimental tank. In particular, the front of the fluid formed by the returning, 
plunging wave is no longer 2D. Apart from this problem there was not enough detailed in- 
formation to determine if a particular method was superior. For example, the experiments 
did not provide them with the decay of the kinetic energy, the decay of the enstrophy, the 
vortex structure or the velocity correlation functions. 

Over the last few years experiments and detailed spectral method simulations have been 
applied to study 2D turbulence in fluids with no-slip boundaries. These include the discovery 
of spontaneous spin-upP, experiments on 2D turbulence in a stratified fluid 9 , and the effects 
of solid boundaries^. These results provide a convenient framework within which to study 
SPH turbulence and this has been exploited by^, ancP to test SPH without turbulence 
models. In addition MonagharP discusses and applies an SPH turbulence model for these 
problems. 

A key feature of turbulence in a square no-slip boundary is that the angular momentum 
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is not constant, and may change abruptly from the initial state then decay (ISHIl). The SPH 
simulations oP predict the same qualitative results though the details depend on the initial 
state. The SPH and spectral methods give results for the decay of the kinetic energy, and 
the enstrophy which are in satisfactory agreement. There are differences between the two 
SPH codes since that of Monaghan and Robinson models the boundary with layers of fixed 
fluid particles, and uses a cubic spline kernel, while that of Monaghan uses boundary force 
particles, and one of the Wendland kernels. Robinson (private communication) found that 
the results for the decay of the kinetic energy converged more rapidly when the cubic spline 
kernel was replaced by the Wendland kernel. However, as we shall show, most of the features 
of the turbulence in a fixed box, for example the kinetic energy decay, the enstrophy decay, 
and the structure of the evolved vorticity field, are very similar. 

In the present paper the study of SPH simulation of turbulence will be extended using 
the SPH code of Monaghan (2011) but without the turbulence. We first confirm that the 
SPH simulation for the case of decaying turbulence in a no-slip square box converges and the 
spin-up is consistent with that found using the spectral theory. Second we study decaying 
turbulence when the box containing the fluid is allowed rotate under the surface stresses 
produced by the fluid. In this case the total angular momentum of the system of box 
and fluid is conserved. Third we simulate turbulence when the fluid and box are in rigid 
rotation when the turbulence is initiated. The evolution of the turbulence was then studied 
both when the box was forced to rotate at its initial angular velocity, and when it was 
allowed to be driven by the fluid stress. This problem is related to turbulence in the earth's 
atmosphere, but a more complete discussion along the lines of the /3-plane study of^, will 
not be attempted. 



II. THE EQUATIONS OF MOTION IN TWO DIMENSIONS 

We consider an incompressible fluid is moving in two dimensions within a square boundary 
with no-slip boundary conditions. It is convenient, especially when we give the boundary a 
mass and allow it to rotate, to refer to the boundary as a box. The acceleration equation 
for the fluid is 

at at P P 
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From this equation it is straightforward to show that the rate of change of the kinetic energy 
E K of the fluid is given by 

where uz is the vorticity, and z is a unit vector perpendicular to the plane of the fluid. The 
integration over a volume is equivalent in the present case to an integration over the two 
dimensions of the fluid, and is denoted by A. Because of the boundary condition, we have 
used the fact that v is zero on the boundary. The rate of change of the vorticity is given by 

- = -,W (3) 

The total enstrophy varies with time according to 

/ \bj 2 dA = -v I (Vu) 2 dA + v I u(Vu ■ n)ds, (4) 

J A 2 J A JB 



d_ 

dt 



where v is the kinematic viscosity, n is an outward unit vector, and the second integration is 
around the boundary. These equations, together with one further equation giving the time 
variation of the gradients of the vorticity, were used by^l to develop his celebrated argument 
for the existence of a turbulent energy spectrum associated with a vorticity cascade to shorter 
length scales. 

The rate of change of total angular momentum L of the two dimensional fluid is given by 



d 

dl 



/ p((r x v) ■ z)dA = / Pr ds + /i / rxV\dA (5) 
J A JB J A 



The second integration can be written in terms of a surface stress, but it is convenient to 
work with the form given noting that for an incompressible fluid 

V 2 v = -V x uz, (6) 

and for a two dimensional fluid 

r x (V x wz) = V(r • wz) — r • Vw. (7) 



Making use of these various relations we find 

dL 
~dt 



Pr ■ ds + / tor ■ nds -2 udA. (8) 

B JB J A 
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The last term can be written as the circulation around the boundary and this vanishes 
because of the no-slip condition. 

In this paper we also consider turbulence in a two dimensional fluid contained within a 
no-slip, square boundary rigidly rotating with angular velocity Qz. In this case it is useful 
to consider the equations of motion in a frame rotating with angular velocity Q. For this 
purpose, we let v denote the velocity in this rotating frame, and wz its vorticity. The 
acceleration equation is then 



where d/dt denotes the derivative following the motion in the rotating frame. By taking 
the curl of this equation we get the same equations for the rate of change of vorticity as for 
the non rotating box. The equations for the rate of change of the total energy and total 
enstrophy, calculated in the rotating frame, are also the same as for the non rotating box. 
We believe it is therefore reasonable to expect that Batchelor's argument would apply when 
the boundary rotates rigidly throughout the simulation. We find that this conjecture is 
correct. 

We also study the case where the boundary is given a mass and a moment of inertia (we 
then call it a box) and allowed to rotate freely under the stresses from the fluid from the 
moment the turbulence is initiated. In this case the total angular momentum of the fluid 
and the box is constant. In addition, the stress felt by the fluid is less because the box moves 
in response to the stress on it from the fluid. The effect is similar to a person attempting 
to walk on a platform that is free to move. One consequence of this is that the change in 
the angular momentum of the fluid is a factor 1/500 less than when the box is fixed. Our 
simulations also show that the form of the vortices that evolve is very different from that 
in the case of a fixed box, or a box which rotates rigidly throughout the simulation. There 
is no simple expression for the rate of change of the fluid energy and vorticity in this case 
because the dynamics of the fluid must include the effects of the motion of the box. 



~dt 



= -2Qz x v + Q 2 r - 




(9) 



III. SPH MODEL 



We consider a weakly compressible fluid with pressure P a function of density p. Surface 
tension is neglected. The reader is assumed to be familiar with standard SPH as described in 
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the reviews by* 14 * 15 l In the following, the labels a and j are used for SPH fluid and boundary 
particles respectively, and rj is used when a summation is over both fluid and boundary 
particles. The SPH form of continuity equation is 

d ^=PaJ2 — (Va-Vr > )-V a W arp (10) 

where the mass, position, velocity, density and pressure of particle a are m a , r a , v a , p a , and 
P a , respectively. The summation is over all particles. The function W ar) = W(\r ar] \, h) is the 
SPH kernel, \r ari \ is the distance between particle a and particle 77, and h = (h a + h v )/2 is 
the average smoothing length. In the calculations to be described the kernel is the fourth- 
order Wendland function^ for two dimensions. This function, when normalized so that 
2n j W(r, h)rdr = 1, is given by 

W(zJ) = ^(2 - z/h)\l + 2z/h), (11) 

if z < 2h, and zero otherwise. Simulations of a wide variety of problems show that the 
choice h v = 1.55, where 5 is the initial particle spacing, gives good results. The interaction 
between any two fluid particles is zero beyond 35. The gradient taken with respect to the 
coordinates of particle a is denoted by V a . The pressure of fluid particle a is given by 

- l) , (12) 




7 \\Po. 

where po is the reference density of the fluid. The speed of sound c s is 10 times the maximum 
speed of fluid V max which we estimate from the initial velocity field. The boundary particles 
have zero pressure. 

The acceleration equation for the SPH particle a is 
dv a 

s 



dt 



„ V ra rrt / „■ 



The first summation in ( 13 ) is over all particles and the second is over the boundary particles. 



The viscosity is determined by U ari for which we use the form 

n av = ~a^^f^, (14) 

Par] I "ar] | 

where a is a constant, v a?? = v a — v^, r ari = r a — r^, p arj = (p a + Pr,)/2 denotes the average 
density, and c = (c a + c v )/2. The constant a can be written in terms of the kinematic 
viscosity by converting the summations to integrals. We find for the Wendland kerneP- 7 that 
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-ahc. (15) 



The last term in Eq. (13) is the boundary force on fluid particle a due to the boundary 
particles. This force, together with the viscous forces due to the boundary particles included 
in the first term, is equivalent to the Siroviclf^ formulation of the boundary conditions in 
terms of boundary forces, and closely related to the Immersed Boundary Method of PeskinP^. 



f a j is given by^ 

f . = a 'J i^L Mr) 

aJ (|r aj |-A)|r a ,|' 

where § a j = ^(l + §<? + 2g 2 )(2 — q) 5 for q < 2 and is otherwise 0. T is a constant equal to 
2V^ ax / \m a + rrij) , A = 5/3 is the boundary particle spacing , and q = \r a j\/h a j. 

The position of any fluid particle a is found by integrating 

t - - < 17 > 

For convenience in describing the time stepping algorithm we write the equations in the 
form 

^ = F a (r,p,v), (18) 

and 

£ = D.(r,v). (19) 

These equations were integrated using a time stepping scheme that is second order and 
based on Verlet symplectic method. In the following A denotes a quantity A at the begin- 
ning of the current time step, A 1 / 2 at the midpoint of the step, and A 1 at the end of the 
step. The time stepping equations, where St is the time step, can then be written 

*l /2 = r° a + f v°, (20) 
vy 2 = *(v° + fF°), (21) 
PT = P° a + f D°. (22) 

where k is a damping factor used to bring the particles to equilibrium before initializing 
the turbulence. This is required because the boundary forces are initially unbalanced. The 
simple and efficient method which we use here is that k is 1.0 for all steps except that every 
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4th time step it is set to 0.9854. This value of k is optimal to reduce the kinetic energy to the 
less than 0.01 percent of the total kinetic energy of vortices after 1000 damping steps. For 
each simulation a minimum time is necessary for damping, so that the number of damping 
steps will increase as resolution increases (5 made smaller). 

The time step is completed by calculating v a , r*, and p a according to 



(v° + 5tFy 2 ), (23) 
rl = + fv a , (24) 

pi = pi /2 + yDi(ri, P ya,vi). (25) 

To improve the speed we replace F° in 22 by Fa 1//2 , i.e half a step back. 
SI units are used throughout this paper. 



IV. THE INITIAL TURBULENT VELOCITY FIELD 



For the numerical study of 2D decaying turbulence in a container the initial velocity field 
has been set using either Chebyshev polynomials^^!! or Gaussian vortices^ES. The results 
for both setups are similar. In this paper Gaussian vortices are used. For the general case 
NxN equal-size Gaussian vortices are placed on a regular lattice, in a checker-board pattern 
of positive and negative vorticity. The initial distance between the centres of vortices, A is 
S/(N + 1), where S = 1 is the width of box. The centres of the vortices were then given a 
random displacement 0.06A(2r — 1), where r is a quasi random number between and 1. 
The rule for the random number is jran = mod(jran * k + I, m) and r = jran/rrW^. Here 
k, I, and m are 106, 1283 and 6075, respectively. When we wish to compare the effect of 
different initial states we use 12 different initial values of jran. The values were 11, 13, 17, 
23, 37, 49, 137, 191, 1111, 1117, 1139, and 3737. 

The SPH particles were placed on a grid of squares in the domain < x < 1 and 
< y < 1. Each SPH fluid particle's velocity was calculated by summation over the vortices 
according td-^ 

NxN 

*• - E -^Sf t 1 - e """ 3/ °") • (26) 

k=l 1 aKl 

where, is the center of the kth vortex, Uk is the vortex strength and its absolute amount 
is equal to 1.0 for all vortices, and is the amplitude of vortex k equal to S/AN. Finally, to 
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ensure that the no-slip condition is closely satisfied in the initial state, a smoothing function 
similar to that used by other authors^ was applied . The final initial velocity field for fluid 
particle a is v a (x,y) = f(x a )f(y a )V a where f(x) = 1 — exp(—9x' 2 ) and x' = 1 — (2x — l) 2 . 

V. FIXED BOX: ENERGY, VORTICITY AND SPIN UP 
A. General results 

In order to confirm that our SPH code behaves correctly it is convenient to compare our 
results against the experimental and computational results $S!MMMS1 They showed that 
the energy and enstrophy decayed approximately as a power of t until the flow is dominated 
by viscosity, and the angular momentum of the fluid may first increase in magnitude then 
slowly decay. This change in angular momentum is due to the reactive stress exerted on 
the fluid through its interaction with the walls. To show that these phenomena are also 
produced by our code we simulated turbulence in a fixed square of side lm. The turbulence 
was initiated using 10 x 10 Gaussian vortices and 275 x 275 SPH particles and = 2000 as 
described in the previous section. The SPH particles were damped for 6000 time steps after 
which the turbulence was initiated. 

The vorticity uj a for particle a, was determined by assuming that velocity of the particles 
relative to particle a was a linear function of the relative coordinates x and y. The coefficients 
of this function were determined by least squares after which the derivatives of the velocity 
components were calculated. In the figures the vorticity field has been scaled between and 
1. The total enstrophy of the flow was then calculated according to 

( 27 ) 

* Pa 
a 

Fig. [I] shows that the average size of the vortices increases with time as found by^ ancPl 
During this time intense shearing occurs near the boundaries as can be seen, for example, 
in the first frame of the second row on the left and right sides. The increase in the absolute 
value of the fluid angular momentum (the spin-up) shown in the Fig. |4j accompanies the 
formation of relatively large vortices. After spin-up, the angular momentum decays slowly, 
and this stage is characterised by a relaxation process to a monopolar structure that is more 
or less situated in the centre of the container. Subsequently, the flow relaxes viscously (the 
last frame in Fig. [T]). 
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Figure 2(a) shows the variation of kinetic energy for the present calculations (using 275 x 
275 particles), those oP, and those oP against scaled time t (time is scaled by S/2Urmsi 
where Urms is the initial RMS velocity of particles and S is the width of box). The scaled 
time for the results oPare based on the best estimate of the time scaling used in their graphs. 
All of these calculations indicate that E% ~ 1/t 8 for < t < 20, with a more rapid decay 
for greater t. The agreement between the three sets of results is very satisfactory. In figure 
2(b) we show the results of the experiments and the present SPH results for the variation 
of the normalized kinetic energy with time. The agreement is again very satisfactory. The 
decay of enstrophy is shown in figure 2(c). The results for the different numerical codes 
and the experimental results are similar to each other, and show an enstrophy decay that 
is faster than the decay of the kinetic energy and given by ~ l/t lA . We can reasonably 
conclude that the decay of energy and enstrophy is not sensitive to the details of the initial 
velocity field (Gaussian vortices or a Chebyshev expansion), or sensitive to the two different 
treatment of boundaries in the SPH simulations. 



B. Convergence 

In Fig. [3^i shows the variation of E k with t for different values of n y the number of particle 
spacings along a side of the boundary. The initial particle spacing is therefore l/n y . The 
graphs for E k indicate that the results are converging. In Fig. [3}} we show E k at times 5, 
10, 15, 20 and 25s as a function of resolution. These results indicate linear convergence. 
However, in order to determine the convergence of any quantity A we follow^, and assume 

A(S) = A(0) + P5 p , (28) 

where A(0) is the value at infinite resolution, 5 is the initial particle spacing, and (5 and p are 
parameters to be calculated. The parameters A(0), (3 and p can be found by evaluating A 
at three values of 5. The trend of the kinetic energy with resolution, which is nearly linear, 
is shown in Fig. [3j 

Fig. [3}b, shows that for all times the change of kinetic energy with particle size, 5, is 



close to linear for high resolutions. From the expression (28) with Ek for n y = 100, 200 
and 250, we can calculate Ek(Q), (3, and p at other times , The resulting values are given 
in Table 1. Using the calculated values of E k (0), (3, and p, it is possible to estimate E k at 
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FIG. 1. Vorticity plots of 2D decaying turbulence flow in a square container with no-slip boundary 
conditions (JR. = 2000). The vorticity has been scaled between and 1 in order to show the field 
at the later times, t = 0,8.08,13.5,27,35.14,48.7,73.1,146.5 and 309.5 respectively from left to right 
and top to bottom. 

other times and other resolutions. Thus in Table I, the two columns under 275x275 show 
results calculated using SPH (the left column), and results using the convergence formula, 
denoted by E est , are shown in the right column. The agreement between the two columns 
is very good. In general the difference between the SPH values, and those estimated from 
the convergence formula, decreases with increased resolution. 

Data presented in this table show that the kinetic energy converges less rapidly as the 
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FIG. 2. (a) The normalized energy, E/E(t = 0) against scaled time for 3ft = 1500. Note that 
the results of the present calculation and those of Robinson and Monaghan (2011) and Clercx 
et al. (1999) are very similar, (b) shows the decay of normalized kinetic energy for the present 
calculations and the experiments of Maassen et. al. (2002) for K = 2000. Note the change in the 
maximum and minimum time. The decay of the normalized enstrophy, = 0) for 3ft = 2000, is 

shown in (c) for the experiments, the present simulations, and those of Clercx et al. 
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FIG. 3. The left frame shows the decay of kinetic energy with time for different resolutions at 
3ft = 1000. The initial particle spacing is given by l/n y . Once n y > 200 the decay curves have 
a similar rate of decay. The right hand frame shows E^ against the initial particle spacing for 
different times. The convergence is close to linear. 
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TABLE I. Comparison of the convergence parameters, the SPH values for kinetic energy together 
with the estimated values Est., from the convergence study at 3? = 1000. The resolution is given 
by the numbers at the head of the 5th and later columns. 



t 


E k (t) 


-/3 


P 


275 x 275 


225 x 225 


150 x 150 


75 x 75 










SPH 


Est. 


SPH 


Est. 


SPH 


Est. 


SPH 


Est. 


5 


1.255 


247.4 


1.475 


1.192 


1.193 


1.179 


1.171 


1.093 


1.103 


0.957 


0.831 


10 


1.05 


36.33 


0.962 


0.884 


0.886 


0.855 


0.851 


0.746 


0.756 


0.524 


0.476 


15 


0.845 


29.96 


0.917 


0.666 


0.671 


0.638 


0.636 


0.527 


0.542 


0.319 


0.272 


20 


0.730 


20.04 


0.831 


0.541 


0.541 


0.506 


0.506 


0.398 


0.417 


0.212 


0.173 


30 


0.632 


6.786 


0.585 


0.379 


0.377 


0.343 


0.346 


0.252 


0.269 


0.111 


0.088 



time increases, and for longer times, the order of the convergence is less than one. The 
reason for the absence of second order convergence is not clear to us. 
These convergence results indicate consistency 



C. Angular momentum 



As remarked earlier, a feature of the turbulence in two dimensions is the increase in the 
magnitude of the angular momentum, the spin-up. To confirm that our SPH simulation 
recovers the spin-up, we used 4x4 vortices with 150 x 150 particles and 9ft = 1000. As 



noted earlier in §iV| we choose 12 different values of jran. The time variation of the angular 
momentum is shown in Fig. [4| Among these 12 simulations, 5 show a rapid, spontaneous 
spin-up, and a further 4 show spontaneous spin-down, which always characterised by the 
presence of a strong monopolar or a rotating tripolar structure. The other 3 show no spin-up 
or very slow spin- up. During the intermediate stage of these runs a dipolar or quadrupolar 
structure is usually found. For example, the velocity field after 206 seconds is shown in 
Fig. [5j for two different initial setups. This left frame of this figure shows the formation 
of a monopolar vortex which corresponds to the sudden increase in angular momentum. 
The right frame shows the formation of dipolar vortices, which corresponds to nearly zero 
angular momentum during the decay of the turbulent flow, ^reported similar behaviour for 



14 




2 t 4 



50 



100 



150 



200 



FIG. 4. The right frame shows the angular momentum against time (in second) for different random 
positions of the vortices forming the initial state in a fixed box at 5ft = 1000. The left frame shows 
a magnified view of initial angular momentum. 




FIG. 5. Velocity field at t = 206s, for two different random initial setups in a fixed box at 
K = 1000. Left frame jran = 2239 shows a strong negative angular momentum, while the right 
frame calculated with jran = 1117 shows negligible change in the angular momentum. 



an initial velocity field of Gaussian vortices or Chebyshev polynomials as dicP for Chebyshev 
polynomials. 

Although the change in angular momentum is considerable when jran is changed the 
corresponding difference in the decay of kinetic energy with time is small, and all runs have 
a similar power law decay. 
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D. Summary of the results with a fixed box 

Our SPH simulation recovers the principal phenomena found by other authors for turbu- 
lence in a two dimensional box with no-slip boundaries. This includes the growth of vortices, 
the decay of the kinetic energy and enstrophy, and the initial rapid increase in magnitude of 
the angular momentum followed by a slow decay. We conclude from this that our SPH code 
gives satisfactory results for two dimensional turbulence in a no-slip box and, in particular, 
the effectiveness of our boundary force model is confirmed. 



VI. FREELY ROTATING BOXES 



We now consider the fluid angular momentum when the turbulence is initiated and the 
box is allowed to rotate under the torque produced by the fluid. In this case the angular 
momentum of the box and fluid is constant. In order to amplify the rotation of the box we 
give it a mass = 0.001 x (mass of fluid) and moment of inertia /box — \ m boxS 2 ■ The 
force fj on boundary particle j due to fluid particles is given by 

»'Ar (29) 

a 

and the torque on the box about its centre is 

T = ^r io xf i , (30) 

3 

where r JO is the coordinate vector of boundary particle b, relative to the centre of the box. 
The angular velocity of box Q is calculated from 

4-§-T. (3!) 

The initial set up is same as the previous section. Fig. [6] shows the angular momentum 
for the different values of jran. The left frame of Fig. [6] shows the angular momentum of 
the fluid. The key point is that the fluid angular momentum remains nearly constant for 
all values of jran. The right hand frame of Fig. [6] shows the angular momentum of the 
box which may be positive or negative but the magnitude is always small and decreases 
with time. The magnitude of the change in the fluid angular momentum can be calculated 
from the box angular momentum. In the case of the fixed box the fluid angular momentum 
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FIG. 6. The effect of initial total angular momentum of fluid particles on a freely rotating box 
at 3? = 1000. a) total angular momentum of fluid particles, b) total angular momentum of 
box particles. Total angular momentum of box plus total angular momentum of fluid particles is 
constant for each set. (the numbers in legend show jrarfs in producing random numbers algorithm) 

typically changes by ~ 1, whereas in the present case it typically varies by 0.002, a factor 
500 less. The average size of the vortices increases with time as in the fixed box. For all 12 
initial setups, a monopolar vortex formed at the centre of box. In contrast, some of the 12 
setups for the fixed box, produce dipole or quadropole vortices, e.g. compare Fig. [7] with 
Fig. |5j The decay of enstrophy in the freely rotating box is similar to that in the fixed box. 

The decay of total kinetic energy of the fluid in the rotating box and the fixed box are very 
close, although, in the last stages of rotation, when the kinetic energy is ~ 10~ 3 , the kinetic 
energy of the fluid in the rotating box is larger than the kinetic energy for the fixed box. 
In Fig. [8] the results of different runs are sorted according to the initial angular momentum. 
For each run the initial and final kinetic energy, final box angular momentum, and final 
box rotation are shown. The kinetic energy is multiplied by 10 4 in this figure and the box 
angular momentum is multiplied by 10 3 . The direction of rotation of the box is given by the 
sign of initial angular momentum of fluid, but the magnitude of box angular velocity is not 
related to the magnitude of the initial angular momentum. The initial kinetic energy for all 
runs is nearly same, but the kinetic energy after 206 seconds can be significantly different 
according to the initial set of random shifts of the vortices. 
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FIG. 7. Velocity field at t = 206s, for two different random initial setups in a freely rotating box 
at K = 1000. Left frame jran = 2239 shows a strong negative angular momentum, while the right 
frame calculated with jran = 1117 shows a strong positive angular momentum. 




FIG. 8. Initial total kinetic energy, Ek(0), final total kinetic energy, total fluid angular momentum, 
L fluid, total box angular momentum, Lf, ox , and total box rotation, 9, after 206s for a free-rotating 
box at K = 1000. The horizontal axis shows the number of initial setups for different random 
numbers, which are sorted according to Lfi ui ^. 



VII. THE EFFECT OF BACKGROUND ROTATION 

In the previous cases the initial angular velocity f^o of box and fluid before inserting the 
Gaussian vortices was zero. Another interesting case that is relevant to turbulence in the 
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earth's atmosphere or ocean, is the behaviour of turbulence when the fluid has a background 
rotation. Forced quasi-two-dimensional turbulence in rotating containers has been studied 
both by experiment and by numerical simulations^^. By adding an oscillation to the 
rotation of the container it is possible to determine how shear near the boundary generates 
vorticity. The vorticity generated in this way subsequently moves into the central regions 
of the container though in the SPH simulations oP, the vorticity generated in this way 
dissipated more rapidly than in the calculations oP^. 

In this section we study the turbulence when the fluid and box are set rigidly rotating 
with angular velocity Qq, after which the vortices are added. The vorticity of the initial state 
is then u = 2Q . The additional vorticity due to the gaussian vortices satisfies the usual 
vorticity diffusion equation for two dimensions (equation [3]). If the turbulence is studied in 
the frame rotating with angular velocity Qq the centrifugal force is balanced by the radial 
pressure gradient. The Coriolis force, 2(Qqz) x v, changes the motion of the flow arising 
from the vortices. 

The previous Gaussian vortices were used with jran = 11, and the rotation was chosen 
so that the Rossby numbers, Ro = 2U rms /floS were in the range 0.57 < R Q < 11.0. U rms is 
the root mean square velocity of vortices, which is equal to 0.057, and S = 1 is the width of 
box. The results for infinite Ro were discussed in the previous section. Here Ro ~ 11, 1.1, 
and 0.57 equivalent to Q = 0.01, 0.1, and 0.2 respectively. The Reynolds number was set 
to 1000 for all studies hereafter. 

After damping, the box and fluid were rotated with f2o for 2000 time steps after which the 
box and fluid were rotating rigidly. The Gaussian vortices were then added to the velocity 
field. The subsequent velocity and vorticity field then depends on the rotation of the box 
and how the added vortices change in response to the background rotation. We consider 
two cases. The first, denoted by A, is where the rigid rotation of the box is maintained. The 
second, denoted by B, is when the box is allowed to respond freely to the fluid stresses on 
it. In the following we denote the kinetic energy and the enstrophy in the inertial frame by 
E k and £ respectively, and in the rotating frame by E k and ^ rot . 
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A. The box with maintained rigid rotation 

In this case the box rotates rigidly with angular velocity Jl . The total kinetic energy of 
box and fluid is shown in the left frame of Fig. [9] for Qq = 0.2. The dashed line shows the 
kinetic energy for the rigidly rotating box. The kinetic energy increases suddenly because 
of the extra kinetic energy of the vorticity field. This extra kinetic energy decays rapidly to 
a steady state where the fluid and box are rotating rigidly. The relative decay of the kinetic 
energy, E k , due to the vortices at different f2 's is calculated by AE k = (E k — E in )/E in x 100, 
where E in is the initial kinetic energy before adding vortices, and E in is the total initial 
kinetic energy of vortices. The kinetic energy of the fluid in the inertial frame is 

E k = E k + -ml + Q Q L, (32) 

where, E k and L are the kinetic energy and angular momentum of the fluid in the rotating 
frame, respectively, and / is the moment of inertia of the fluid around the box centre. E k 
and L are zero before adding vortices. 

The calculated kinetic energy in the rotating frame is the same for all JTs, for a given 
set of vortices. The values of AE k , using E k as the minimum kinetic energy which occurs 
nearly at time 100s, is shown in column 7 of table [TTJ The results in table [TT] show that the 
rate of decay of the kinetic energy of the fluid increases significantly by increasing the box 
angular velocity. Therefore, while the rate of change of kinetic energy in the rotating frame 
is always negative, the rate of change of kinetic energy in the inertial frame can be negative 



or positive, due to the last term in equation 32 This can be seen easily from the left frame 
of Fig. [9} 

The right frame of Fig. [9] shows the total angular momentum of the box and the fluid for 
the case where the box rotates rigidly. The angular momentum of fluid in the inertial frame 
decreases to a minimum (column 3 of Tab. [TT]), then it decays gradually until at infinity it 
reaches the initial angular momentum of fluid before adding Gaussian vortices (column 4 of 
Tab. [TT]). Recall that in this section the box rotates rigidly throughout the simulation. The 
change in the angular momentum of the fluid is almost independent of Qq, e.g the maximum 
change in the angular momentum of fluid is shown in column 5 of Tab. [TTJ which shows the 



subtraction of column 3 from column 4. In the left frame of Fig. [10] we show the decay of 
kinetic energy in the rotating frame for different values of the angular velocity of the box. 
As conjectured after (9) it is expected that the decay of the kinetic energy calculated in the 
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TABLE II. The rate of decay in kinetic energy, AE k , total angular momentum, of box Lbox-, an d 
fluid, L Fluid for when the box continues to rotate rigidly at different angular velocities, Orj. Ei n is 
the initial kinetic energy of vortices. 



LFluid 



(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


Q 


Lbox 


Min 


Max 


AL 


E- 


AE k (t = 100) 


0.00 


0.000 


-0.94 


0.00 


0.94 


1.575 


+0.092 


0.01 


0.002 


0.414 


1.61 


1.20 


1.575 


-0.202 


0.10 


0.017 


14.84 


16.10 


1.26 


1.573 


-6.331 


0.20 


0.033 


30.87 


32.21 


1.35 


1.575 


-14.63 
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FIG. 9. Total kinetic energy, Ek (left) and total angular momentum of box, Lb and fluid, hp 
(right) for different background angular velocities, f2o- "forced" denotes the case of box with 
maintained rigid rotation, and "free" indicates that the box is allowed to rotate freely after the 
vortices are initiated. In the right hand frame the dashed and thin continuous lines denote the 
angular momentum of the fluid for the forced and freely rotated box respectively, and the dotted 
and thick solid lines denote the angular momentum of the box for the forced and freely rotated 
box respectively, calculated in the inertial frame. 



rotating frame should be independent of the rotation because the equations for the decay 



are identical to those in a fixed box. The right hand frame of Fig. [10] shows that the decay 
of Ek varies significantly with the angular velocity of the box. 

The enstrophy calculated in the rotating frame, £ ro4 , is independent of VLq and is same 
as the fixed box, this is shown in the right frame of Fig. [TTJ Using the no-slip boundary 
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Time (second) Time ( second) 



FIG. 10. The left frame shows the decay of the kinetic energy in the rotating frame for a rigidly 
rotating box. As conjectured it is independent of the angular velocity. The right frame shows the 
change in the kinetic energy with time as measured in the inertial frame. 



conditions it can be shown that the enstrophy £ calculated in the inertial frame is 8QlS 2 
larger than £ rot , where S is the half-width of box. As a consequence, when the enstrophy in 



the rotating frame has decayed to zero, £ = 8QqS 2 as shown in the left frame of Fig. 11 graphs 
labelled 'forced '). 



B. The box freely-rotating after the vortices are initiated 

In this case the box rotates freely after the vortices are initiated in the rigidly rotating 
fluid. The subsequent motion of the box is due to the stresses on it from the fluid. This class 
of problem is relevant to a wide range of engineering fluid-structure interaction problems 
such as vortex-induced vibrations. The total kinetic energy of box and fluid in the inertial 
frame, Ek, for this case is shown in the left frame of Fig. [9j By adding Gaussian vortices 
to the background velocity field, a sudden rise in kinetic energy occurs which is followed by 
a sharp decline to the background kinetic energy. The kinetic energy is calculated at the 
same times that the kinetic energy of the box with maintained rigid rotation is minimum 



for each Q,q. The relative decay is shown in the last column of table III The decay rate of 
total kinetic energy (liquid and box) when the box is freely rotating is less than the when 
it is rigidly rotated, and this difference increases significantly by increasing the box angular 



velocity (see Tab. [XT] and Tab. Ill), around 15% for £7 = 0.2. The decay rate of the kinetic 
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FIG. 11. The effect of angular velocity, Uq, on enstrophy, £, a) total enstrophy, calculated in 
inertial frame, b) net enstrophy due to Gaussian vortices, calculated in rotating frame, "forced" 
denotes the case of box with maintained rigid rotation, and "free" indicates that the box is allowed 
to rotate freely after the vortices are initiated. £o is the enstrophy of initial Gaussian vortices equal 
to 0.999. 



TABLE III. The rate of decay of kinetic energy, AE^, total angular momentum, of box Lbox, and 
fluid, L Fi u id for the case of a box which is rotated initially with Qq then rotates freely after the 
vortices are initiated. 







LBox 






LFluid 






(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


^0 


Min 


Max 


AL 


Min 


Max 


AL 


AE k (t = 100) 


0.00 


0.0 


0.005 


0.0054 


-0.126 


0.0 


0.126 


+1.84 


0.01 


0.002 


0.007 


0.0051 


1.484 


1.609 


0.126 


+1.67 


0.10 


0.017 


0.022 


0.0053 


15.963 


16.107 


0.144 


+1.46 


0.20 


0.033 


0.038 


0.0050 


32.060 


32.268 


0.208 


+0.22 



energy increases as the angular velocity increases for both the freely rotating box and the 
rigidly rotating box. 
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The right frame of Fig. [9] shows the total angular momentum of box and the total angular 
momentum of fluid, when the box is free to rotate after the gaussian vortices are initiated. 
The angular momentum of the box is not constant and increases to a maximum, then it 
decays gradually until it reaches the angular momentum of the box before initiating the 
Gaussian vortices. Note that the combined angular momentum of fluid and box is constant 
and the angular momentum of the box is scaled up by a factor 10 3 . The increase in total 



angular momentum of the box is not related to the f2 (column 4 of Tab. III). The initial 
angular momentum of vortices is around -0.126. The thin solid line in the right frame of 
Fig. [9j which shows the total angular momentum of the fluid in the inertial frame, shows 
a small drop by -0.126 shortly after t — 1, then remains nearly constant, varying slightly 
as the box angular momentum changes. These changes in the angular momentum are in- 



dicated in columns 5 to 7 of Tab. |III| Thence the angular momentum of the fluid remains 
approximately constant for each Q . 



From the right frame of Fig. 11, and focussing on the graphs labelled 'free', it can be 
seen that the decay of enstrophy is independent of Qq as in the case of the rigidly rotating 
box. The enstrophy of the fluid when the box is freely rotated is always slightly larger than 



when it is rotating rigidly. The right frame of Fig. 11 shows the decay of the enstrophy after 
subtracting the initial background rotation from the total flow. Although the rotation of 
the box can change when it is freely rotating it is convenient to consider the enstrophy in 
the frame rotating with the initial angular velocity. The enstrophy shown in the right frame 



of Fig. 11 is calculated in this frame both for the forced and free rotations. Our comparison 
shows that the enstrophy decay in this frame is nearly same for all fVs, and is very close 
(0.17 %) to that of the rigidly rotated box . 



The angular velocity of the box due to the vortices is shown in the left frame of Fig. 12 



This figure shows that the effect of vortices is nearly the same for all background rotations 
and, as a consequence, the total rotation of the box is similar for all values of Q as shown 



in Fig. 12 for different fVs. 
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FIG. 12. The effect of background rotation, Qq, on the box rotation. The box is free to rotate after 
the gaussian vortices are initiated; a) The box angular velocity due to vortices, O — Qq, versus time, 
b) The net rotation angle due to the vortices versus time, which is calculated after subtracting the 
forced rotation angle, Clo.t, from the free rotation angle. 



C. Summary of the results with Background rotation 

In the case of the box with maintained rigid rotation, if the flow is considered in a rotating 
frame, the effect of background rotation on the flow characteristics like vorticity structure 



(see Fig. 13), enstrophy, kinetic energy of vortices, and angular momentum of the fluid is 
negligible. These results are consistent with the fact that the equations for the decay of 
enstrophy and kinetic energy in the rotating frame are identical to those for a fixed box. 
The decay rate of the kinetic energy in the inertial frame increases with Qq. 

When the box is allowed to rotate freely, the effect of the background rotation is negligible, 
though the vortex structure is completely different from that for the box with maintained 



rigid rotation (see Fig. 13). 



VIII. CONCLUSION 

We have described the results of numerical simulations of decaying two-dimensional tur- 
bulence inside a rotating square container with rigid no-slip boundaries, using SPH. Three 
cases were studied here and, in all cases, the turbulence was initiated by superimposing 
the velocity field of a set of gaussian vortices. The first two cases involve a fixed box, and 
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FIG. 13. Vorticity and velocity field for different angular velocities, first, second and third columns 
from the left show results at time=76s after adding vortices for H,q = 0.0, 0.1, and 0.2, respectively. 
The top row show the box with maintained rigid rotation, and the bottom row show the freely 
rotating box results. All velocity fields are calculated the in rotating frame. 

a box that is free to rotate under the stresses from the fluid. The third involves adding 
gaussian vortices to a box and fluid rotating rigidly. After which the box was either forced 
to maintain its rigid rotation, or allowed to rotate under the stresses from the fluid. 

Our simulations confirm that our SPH code correctly simulates turbulent flows in a 
fixed square box with no-slip boundaries. In particular, the simulations show a rapid self- 
organisation of the flow toward one or two large vortices as the total kinetic energy of the 
flow decreases. By changing the random shift of the initial vortices from a regular grid 
the initial conditions can be changed. The resulting simulations provide an ensemble from 
which more general properties can be deduced. It was observed that 5 of them show a rapid 
spontaneous spin-up, and a further 4 showed spontaneous spin-down. All of them were 
accompanied by an increase in the magnitude of the angular momentum of the fluid. That 
behaviour was always accompanied by a strong monopolar or a rotating tripolar structure. 
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The other 3 simulations showed no spin-up or very slow spin-up. During the intermediate 
stage of these simulations a dipolar or quadrupolar structure was usually found, and the 
net angular momentum of flow remains approximately constant or decayed only very slowly. 
The decay of the total kinetic energy of the flow was more rapid in all cases. These results 
are in good agreement with other experimental and numerical studies 2 - 9 2S . 

When the box was allowed to rotate under the stresses from the fluid the angular mo- 
mentum of the fluid remained approximately constant and the total angular momentum of 
box and fluid remained constant as it should. Simulations with slightly different initial con- 
ditions change the dynamics of the box. In particular, the box angular momentum changes 
and therefore the amount and direction of its rotation. All 12 initially different setups con- 
cluded in a strong monopolar vortex. The rate of decay of the kinetic energy is nearly the 
same as that for the box with maintained rigid rotation. 

From the study of the effect of a background rotation on the decaying turbulence, it was 
seen that for the box with maintained rigid rotation, the effect of box angular velocity on 
the angular momentum, vorticity structure, velocity field, kinetic energy, and enstrophy is 
negligible, when these quantities are calculated in a frame rotating with Q Q . This is consistent 
with the fact, already noted, that the equations for the decay of kinetic energy and enstrophy 
in the rotating frame are the same as in the inertial frame. For the box freely-rotating after 
the turbulence was initiated, kinetic energy, angular momentum,vorticity structure, velocity 
field, enstrophy, and rotation, calculated in the rotating frame were independent of flo. 
Although some parameters like kinetic energy or enstrophy calculated in the rotating frame 
are similar for both the freely rotated box and the box with maintained rigid rotation, the 
vorticity structure, velocity field and angular momentum are completely different for these 
two situations. 

The resolution required for these SPH simulations depends on the accuracy required and 
the time for which that accuracy should be maintained. The results of this paper and those 
oP, ancP, show that for the square box of half width S the number of particles should 
be ~ 200 with particle spacing dp = 2S"/200, in order to determine the energy decay to 
within 10% for t < 15. The resolution length can by estimated from h = 1.5dp and this 
is ~ 0.55/v^R which agrees with the estimate oM We also note that our treatment of 
the boundaries using boundary forces gives results in good agreement with those oP for the 
energy decay obtained using a boundary modelled by layers of fixed fluid particles. However, 
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the convergence of our results for the energy decay, obtained using the Wendland kernel, 
appear to be faster than those oP (see for example their figure 9) using the cubic spline 
kernel. This agrees with the experience of Robinson (private communication) who repeated 
some of his calculations using the Wendland kernel instead of the cubic spline kernel. 

In addition to the class of problems considered here SPH can be applied to other turbu- 
lent flows where it has many advantages. One of these is the turbulent flow produced by 
physical stirrers, for example cylindrical rods moving on specified paths. Such a problem 
can be simulated easily with SPH and, in unpublished work, we have studied the turbulence 
produced by such stirrers moving on a variety of trajectories. SPH has also proven useful for 
studying breaking waves especially those formed by a flow hitting and running up a wall and 
finally forming a backward breaking wave. A region of strong turbulence is created by the 
impact of the breaking wave on the incoming fluid. This problem, which is a key feature in 
sloshing in marine tanks, could be tackled using an SPH code of the kind we have described. 
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